#include "field.hpp"
#include "meshTypeFunctions.hpp"
#include "bc_orthogonalFlowVelocity_velocity.hpp"
#include "bc_orthogonalFlowVelocity_waterdepth.hpp"
#include "util.hpp"

#include "testHelpers.hpp"

#include <catch2/catch_test_macros.hpp>
#include <catch2/catch_approx.hpp>

using namespace hms;
using Type = hms::boundaryCondition::Type;
using Catch::Approx;

TEST_CASE( "OrthogonalFlowVelocity_velocity" ){
	auto test = []( auto&& dummyMesh ){
		using MeshType = remove_qualifiers<decltype(dummyMesh)>;
		auto mesh { uniMeshAs<MeshType>(4, 4, 1, 1, true) };
		using BC = OrthogonalFlowVelocity_velocity<MeshType>;

		const BoundaryPatch& left   { mesh.boundaryPatch("left") };
		const BoundaryPatch& right  { mesh.boundaryPatch("right") };
		const BoundaryPatch& bottom { mesh.boundaryPatch("bottom") };
		const BoundaryPatch& top    { mesh.boundaryPatch("top") };

		Field v { mesh, 2 };
		Field d { mesh, 1 };
		scalar vBoundary { 2.0 };

		v.values() = 1.0;
		d.values() = 1.0;

		for ( const BoundaryPatch& p : mesh.boundaryPatches() ){
			v.template setBoundaryCondition<BC>(p, vBoundary).update();
		}

		REQUIRE(
			v.boundaryCondition(left).type()
				== Type::orthogonalFlowVelocity_velocity
		);

		REQUIRE( (v.ghostValues(left).row(0) == 2).all() );
		REQUIRE( (v.ghostValues(left).row(1) == 0).all() );
		
		REQUIRE( (v.ghostValues(right).row(0) == -2).all() );
		REQUIRE( (v.ghostValues(right).row(1) ==  0).all() );

		REQUIRE( (v.ghostValues(bottom).row(0) == 0).all() );
		REQUIRE( (v.ghostValues(bottom).row(1) == 2).all() );

		REQUIRE( (v.ghostValues(top).row(0) ==  0).all() );
		REQUIRE( (v.ghostValues(top).row(1) == -2).all() );
	};
	TEST_ALL_MESHTYPES(test)
}

TEST_CASE( "OrthogonalFlowVelocity_waterdepth" ){
	auto test = []( auto&& dummyMesh ){
		using MeshType = remove_qualifiers<decltype(dummyMesh)>;
		auto mesh { uniMeshAs<MeshType>(4, 4, 1, 1, true) };
		using BC = OrthogonalFlowVelocity_waterdepth<MeshType>;

		const BoundaryPatch& left   { mesh.boundaryPatch("left") };
		const BoundaryPatch& right  { mesh.boundaryPatch("right") };
		const BoundaryPatch& bottom { mesh.boundaryPatch("bottom") };
		const BoundaryPatch& top    { mesh.boundaryPatch("top") };

		Field v { mesh, 2 };
		Field d { mesh, 1 };

		v.values() = 1.0;
		d.values() = 1.0;

		scalar vBoundary { 1.1 };

		for ( const BoundaryPatch& p : mesh.boundaryPatches() ){
			v.ghostValues(p) = vBoundary;
			d.template setBoundaryCondition<BC>(p, v, vBoundary).update();
		}

		scalar orthoVal { std::pow((1-1.1+2*sqrt(9.81*1)), 2.0) / (4*9.821) };

		REQUIRE ( (d.ghostValues(left)   = orthoVal).all() );
		REQUIRE ( (d.ghostValues(right)  = orthoVal).all() );
		REQUIRE ( (d.ghostValues(top)    = orthoVal).all() );
		REQUIRE ( (d.ghostValues(bottom) = orthoVal).all() );
	};
	TEST_ALL_MESHTYPES(test)
}